——— Introducción al análisis espacial ———

¶ Serie de talleres en “Herramientas de análisis cuantitativo y su aplicación en la conservación de la biodiversidad” ¶

Ecodiversa Tropical

1. Objetivos de aprendizaje

Al final de este taller usted podrá:

  • Identificar las características de los datos espaciales, en sus dos formatos: ráster y vectorial.
  • Conocer los paquetes y funciones en R más utilizados para el manejo de los dos formatos.
  • Utilizar bases de datos de libre acceso disponibles en la web.
  • Integrar y graficar los dos formatos
  • Implemetar algunos análisis espaciales que pueden servir de herramientas para la toma de decisiones en la conservación de la biodiversidad.

2. ¿Qué hacer si no tienes conocimiento previo en R?

Si no tiene concocimientos de R le sugerimos tomar estos cursos que le ayudarán a tener bases y desarrollar sus habilidades.

3. Para empezar

Inicie RStudio y abra un nuevo R Script utilizando el menu: File > New File > R Script. Usted podrá copiar los códigos que aquí presentamos y probarlos localmente en su computador. Además habrá unos ejercicios que le permitirán consolidar y aumentar sus habilidades en desarrollar códigos y en el análisis espacial.

4. Paquetes en R que vamos a utilizar

  • terra: Análisis espacial de datos (Hijmans et al 2022).

  • raster: Análisis y modelamiento de datos geográficos (Hijmans et al 2022).

  • sf: Simple Features (en inglés) permite el manejo de datos vectoriales (Pebesma et al 2022).

  • stars: permite el manejo de colecciones espacio-temporales y cubos de datos en formato vectorial y raster (Pebesma et al 2022).

  • gdalcubes: permite el manejo de cubos de datos, por ejemplo observaciones de la tierra a través de colecciones de imágenes de satélite (Appel et al 2022).

  • geodata Paquete para decargar datos geogáficos (Hijmans et al 2022).

raster y terra son paquetes en R que permiten el manejo de datos en formato ráster. terra ha sido desarrollado para reemplazar raster, es mucho más rápido y posee mas funcionalidas (entre ellas el manejo de datos vectoriales). En este vínculo puede ver las principales funciones de terra y la comparación de las principales funciones que se encontraban en raster y cómo se llaman ahora en terra.

Por otro lado, el paquete sf (simple features, en inglés) ha reemplazado completamente el paquete sp. Vamos entonces a enfocarnos en sf que permite el manejo de datos vectoriales, como por ejemplo: puntos, líneas y polígonos.

Al final veremos algunas las aplicaciones de los paquetes stars y gdalcubes que permiten el manejo de cubos de datos; como por ejemplo imágenes de satélite que es una colección de varias bandas y años.

5. ¿Por qué es importante sabe manejar y analizar información espacial?

  • Entender patrones de biodiversidad
  • Predicciones
  • Monitoreo
  • Planeación de la conservación

6. ¿Cuáles son los formatos en que encontramos información espacial?

Los principales formatos en el cual encontramos información espacial son: vectorial y ráster.

Los datos vectoriales pueden representar puntos, líneas o polígonos. En general, este formato es común para almacenar y representar observaciones (ej., presencia de especies), objetos lineales (ej., como vías) y áreas geográficas (ej., límites administrativos).

Los datos en ráster (que se concocen también como malla o grilla) almacenan pixeles (celdas). Imagine una fotografia en su cámara digital, o las imágenes de su televisor, o las obtenidas por los satélites, etc. Estas imágenes están compuestas de pixeles. En adelante utilizaremos el término técnico ráster para hacer referencia a este formato de datos, pero paras dar fluidez al texto tambien puede encontrar malla o grilla. Utilizaremos también las palabras pixel o celda para hacer referencia a la unidad mínima de un ráster que contiene información.

7. Ráster

Como mencionamos, el ráster es un formato de datos que almacena imágenes digitales en forma de pixeles (o celdas). En general es una cuadricula con varias celdas que contienen valores.

Empecemos con unos ejemplos para visualizar y entender mejor de lo que estamos hablando.

Creando rásters

Vamos a crear un objecto ráster de 10 columnas y 10 filas, utilizando la funcion rast() del paquete terra.

r1 <- rast(ncol=10, 
           nrow=10
           )
r1
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84

Al explorar el objeto r1 vemos que:

  • el parámetro class: es SpatRaster (es decir, es básicamente un ráster).
  • el parámentro dimensions: es de 10 filas por 10 columnas (tal y como la definimos) y una capa (nlyr).

Pero además vemos que la función le ha asignado otros atributos como:

  • cood. ref.: la función asigna por defecto coordenadas geograficas lon/lat, usando WGS 84 (World Geodetic System, por sus siglas en inglés).

Ahora bien, la función ha generado únicamente la estructura del ráster, pero todavía el ráster no tiene valores.

Asignemos valores a cada pixel con la funcion values(). Note que la función ncell() calcula el número de pixeles en el ráster. En este caso vamos asignar valores secuenciales.

values(r1) <- 1:ncell(r1)
r1
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     1 
## max value   :   100

La función values() ha asignado varias cosas:

  • los valores de 1 a 100 (min value: y max value:, respectivamente)
  • extent: por defecto la extensión es en coordenadas geográficas. En este caso longitud (-180 a 180) y latitud (-90 a 90; de polo sur a polo norte).
  • resolution: es la resolución de la grilla, es decir el tamaño de cada pixel o celda, se calcula con base en el número de filas y columnas que indicamos al principio (dimensions:) y en el extent:. Es decir, como definimos 10 por 10, la función asigna el siguiente tamaño de pixel 36 x 18 (grados de longitud(x) y latitud(y), respectivamente).
  • el nombre del ráster (name: lyr.1). Lo podemos cambiar más adelante con un nombre que refleje lo que contiene el ráster.
  • el parámetro source: indica la fuente del ráster. En este caso el objeto está en la memoria virtual del computador. Más adelante podrá ver que si descarga un ráster desde su computador podrá utilizar la función sources() para conocer la ruta donde está guardado el ráster.

Visualicemos el ráster que hemos creado, utlizando la función plot() y la función text() para agregar los valores de cada pixel sobre el ráster.

plot(r1) # uilizamos la función plot para visualizar r1
text(r1) # utilizamos la función text para adicionar los valores de cada celda

La grilla que se genera posee valores enteros y la función plot() genera la escala de la barra vertical.

Ya hemos creado nustro primer ráster!

Utilidad de crear un ráster ficticio

El ráster qque hemos creado es lo que conocemos como un raster ficticio (dummy raster en inglés). Las utilidades de construir un ráster ficticio son varias, por ejemplo:

  • podemos utilizar el raster ficticio para probar funciones/operaciones y verificar los resultados antes de aplicarlas a un ráster grande y con muchos datos. Es decir, podemos verificar que nuestras funciones arrojen los resultados esperados, antes de aplicar dichas funciones a bases de datos reales.

  • podemos compatir en la web nuestros códigos para resolver un problema al cual no le encontramos solución. Es decir compartimos un ejemplo reproducible de lo que queremos hacer y así otras personas pueden sugerir modificaciones. Unas lecturas de la importancia de crear un buen ejemplo reproducible: Elio Campiteli, stackoverplow.com).

Me voy a explicar con un ejemplo para el primer caso.

Primero creemos otro ráster, por ejemplo: r2 será igual a r1+1.

Revise el gráfico de abajo. ¿Es lo esperado?

r2 <- r1+1 # r2 es el resultado de sumarle 1 a r1 creada anteriormente

par(mfrow=c(1,2)) #definamos una ventana de una fila y dos columnas para presentar nuestros rásters

plot(r1, legend=F, main="r1")
text(r1)

plot(r2, legend=F, main="r2")
text(r2)

Creemos ahora otro ráster r3 con base en r1 y r2.

La formula que vamos a utilizar es (r1+r2)/10. Es decir, sumar los valores de r1 y r2 y luego dividirlos entre 10.

En el gráfico de abajo se encuentran los resultados.

¿Es lo esperado?

Parece que sí. Por ejemplo, la primera celda es 0.3, que es el resultado de (1+2)/10 = 0.3

Este es un ejemplo muy simple, pero espero que muestre la utilidad de crear rásters ficticios.

r3 <- (r1+r2)/10

plot(r3, legend=F)
text(r3, digits=2)

Apilando rásters

Una de las funcionalidades más interesantes con rásters es que las podemos apilar, es decir crear un objeto con un conjunto de capas .

Por ejemplo, apilemos r1, r2, r3 en un objeto que llamaremos r_stack, utilizando la funcion c() del paquete terra (la funcion análoga en el paquete raster es stack() o brick()).

r_stack <- c(r1,r2,r3)
r_stack
## class       : SpatRaster 
## dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## names       : lyr.1, lyr.1, lyr.1 
## min values  :     1,     2,   0.3 
## max values  :   100,   101,  20.1

Vemos que tenemos ahora un objeto con tres capas que se pueden identificar por el nombre (name:en este caso los nombres son iguales) con sus respectivos valores: min values y max values:.

Dentro de este objeto podemos visualizar cada capa de varias formas:

  • La pila de rásters es como una lista de capas, entonces usted puede acceder a cada una utilizando el nombre del objeto seguido de la capa que quiere entre dos parentesis cuadrados, por ejemplo la capa 2 r_stack[[2]].
r_stack[[2]]
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     2 
## max value   :   101
  • La otra forma es utilizando el signo pesos $ despues del objeto (r_stack$) y le aparecera la lista de las capas que se encuentran en el objeto, identificadas con el nombre (name:). En este caso selecionemos la primera capa (r_stack$lyr.1). El problema es que todas las capas tienen el mismo nombre, entonces usted obtendrá este error:
r_stack$lyr.1
## Error: [subset] you cannot select a layer with a name that is not unique

Aprovechemos y cambiemos los nombres de las capas en esta pila.

names(r_stack) <- c("l1", "l2", "l3") # concadenamos (c) el nombre de cada capa utilizando comillas.

r_stack
## class       : SpatRaster 
## dimensions  : 10, 10, 3  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## names       :  l1,  l2,   l3 
## min values  :   1,   2,  0.3 
## max values  : 100, 101, 20.1

Ahora sí intentemos utilizando el signo $. Parece que si funciona.

r_stack$l2
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        :  l2 
## min value   :   2 
## max value   : 101

Resolución espacial

La resolución espacial de un ráster es básicamente el tamaño del pixel.

Creemos otro ráster r4 de 100 columnas y 100 filas y adicionemos el nombre “l4”.

r4 <- rast(ncol=100, 
           nrow=100,
           )
values(r4) <- 1:ncell(r4)

names(r4) <- "l4" # Note que el objeto es r4  y el nombre de la capa lo cambiamos a "l4"

Comparemos r4 con r_stack$l1.

r4
## class       : SpatRaster 
## dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
## resolution  : 3.6, 1.8  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        :    l4 
## min value   :     1 
## max value   : 10000
r_stack$l1 
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        :  l1 
## min value   :   1 
## max value   : 100

La grilla r4 es mucho mas fina que r_stack$l1 porque:

  • la resolución de r4 es de 3.6 por 1.8. Es decir que cada pixel es un rectangulo con un lado de 3.6 y el otro 1.8.

  • y en el raster r_stack$l1 los pixeles son de mayor tamaño (rectangulares tambien) de 36 por 18.

Note que las celdas no necesariamente tienen que ser cuadradas, especialmente cuando tenemos coordenadas geográficas. Más adelante veremos como los rásters projectados sí tienen celdas cuadradas (bueno, en lo posible). Es decir x, y en resolution: deberían tener valores iguales.

Volviendo al ejemplo anterior, podemos también decir que la resolución de r_stack$l1 es mucho más gruesa. de hecho 10 veces más gruesa que r4.

Como la resolución de r4 es más fina, entonces hay mas celdas (10,000; dimensions: 100 por 100) que en r_stack$l1 (100; dimensions: 10 por 10).

Note que el área seleccionada (extent:) y las coordenadas (coord. ref.:) son las mismas. Visualicemos las dos rásters para ver sus diferencias:

par(mfrow=c(1,2))
plot(r4, legend=F, main="r4, resolución fina")
plot(r_stack$l1, legend=F, main="r_stack$l1, resolución gruesa")

Ahora bien, ¿Podríamos crear una pila de rásters con área (extent) y coordenadas (coord. ref.) similares, pero resoluciones y dimensiones diferentes?

Veamos:

rstack_new <- c(r_stack$l1, r4)
## Error: [rast] number of rows and/or columns do not match

Pues no, el Error nos dice que el número de columnas y filas (dimensions) no coinciden.

¿Qué más podemos hacer con rásters?

Ya vimos algunas cosas que podemos hacer con rásters, como por ejemplo apilarlas. Pero hay algunos métodos importantes cuando se quieren hacer análisis utilizando rásters. Cuando nos referimos a métodos hacemos referencia a una serie de funciones que se pueden aplicar a uno o varios rásters para obtener información relevante en su investigación.

No podemos en este taller ver todas las funciones, pero al final usted tendrá la capacidad de explorar y aplicar más funciones.

Por ejemplo, ya vimos cómo podemos crear un raster utilizando la función rast(), apilar varios rásters con c(), y cómo seleccionar algunas capas dentro de una pila de rásters utilizando [[]] o $.

Ahora veamos algunos métodos y funciones que puden ser mucha utilidad en el análisis espacial.

Métodos locales

Los métodos locales hacen referencia a hacer cálculos entre celdas que se sobreponen en una pila de rásters. Es decir la localización espacial de las celdas es la misma. La figura muestra cómo podríamos calcular la media de dos raster que se sobreponen.

Hay varias funciones para hacer cálculos de este tipo. Por ahora utilizaremos la función app() para demostrar este método.

Veamos el código en R para calcular la media de la pila de tres capas que hemos creado; es decir la funcion app() en terra toma los tres valores de cada pixel de cada capa que se sobreponen y calcula la media (mean). Esto se hace en toda la extensión del área que definimos.

r_stack_media <- app(x = r_stack, # Este es el objeto que representa la pila de rásters
                     fun = mean   # Esta es la función para calcular la media 
                     )
r_stack_media
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 36, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : mean 
## min value   :  1.1 
## max value   : 73.7

Note que el nuevo objeto r_stack_media representa la media de las tres capas en r_stack. R asigna el nombre de la función a la capa final en esta caso name: es mean.

Visualicemos el promedio de esto tres rásters.

plot(r_stack_media)
text(r_stack_media, digits=1, cex=0.8)

Tips!

Note que los decimales están indicados con un punto ., para los miles se utiliza la coma ,. Los países del norte utilizan esta nomenclatura y R tambien.

En Colombia es al contrario (no sabría decir si en toda latinoamérica es igual). Por qué sucede esto? No lo sé.

Lo que sé es que usted tendrá que ajustar las bases de datos para que tengan el mismo formato requerido en R. En especial tablas en excel que usted ha construido y que le gustaría incorporarlas a un análisis espacial.

Usted ya puede ver las potenciales aplicaciones de estas funciones en un pila de rásters. Por ejemplo, calcular la media o la varianza de la temperatura durante un año (una pila de 12 capas) en su área de estudio.

Ahora veamos un ejemplo utilizando datos globales reales!!!!

Datos globales y algunas aplicaciones

Vamos ahora a manejar algunos bases de datos de acceso gratuito (open access, en inglés) que pueden ser utilizadas en análisis espaciales.

Por ejemplo, los datos globales de clima, como temperatura y precipitación, pueden ser utilizados ampliamente para entender los efectos del cambio climático en la distribución de las especies.

WorldClim y CHELSA son dos iniciativas que ofrecen datos de buena calidad a nivel global.

Existen además algunos paquetes en R que permiten descargar esta bases de datos. Por ejemplo, en el paquete raster existe la función getData() que permite descargar datos de WorldClim. Sin embargo, esta función será reemplazada por el paquete geoData. Veamos un ejemplo

Temperatura media anual a nivel global (WorldClim)

Descarguemos la temperatura media anual a nivel global a una resolución gruesa (10 minutos que es ~340 km2) desde WorldClim. Tenga en cuenta que descargará un archivo comprimido (.zip) a su computador.

Existe tambien la posibilidad de descargar las 19 variables disponibles de bioclim utilizando ‘bio’ como variable.

tmean <- geodata::worldclim_global(
    var="tavg", # temperatura media
    res=10,     # resolución gruesa
    path = "C:/Talleres_Ecodiversa_R/Intro_espacial" # folder donde desea guardar los datos
                      )

Puede ver que son 12 capas (una por cada mes, ver el atributo nlyr en dimensions:).

  tmean
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_tavg_01.tif  
##               wc2.1_10m_tavg_02.tif  
##               wc2.1_10m_tavg_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :    -45.8840,   -44.80000,   -57.92575,   -64.19250,   -64.81150,   -64.35825, ... 
## max values  :     34.0095,    32.82425,    32.90950,    34.19375,    36.25325,    38.35550, ...

Visualicemos dos capas (Enero y Agosto)

par(mfrow=c(1,2)) 

terra::plot(tmean$wc2.1_10m_tavg_01, main="Temperatura media - Enero", range=c(-50, 40), mar=c(4, 3, 4, 3), plg=list(shrink=0.9, cex=.8), 
    pax=list(side=1:2, cex.axis=.6))
terra::plot(tmean$wc2.1_10m_tavg_10, main="Temperatura media - Agosto", range=c(-50, 40), mar=c(4, 2, 4, 4), legend=FALSE, 
    pax=list(side=c(1,4), cex.axis=.6))

Excelente, ya vemos algunos rásrters con datos reales. Ahora tratemos de aplicar algunas funciones.

Algunas funciones locales

Vamos a aplicar las funciones stdev()y app() al objeto tmean que tiene 12 capas.

  • Medidas de dispersion stdev()

Veamos ahora en donde hay mayor variación de la temperatura durante el año, utilizando la función stdv() del paquete terra.

tmean_sd <- stdev(x = tmean, # SpatRaster object
                  pop = T, # Si es verdadero (T), entonces la desviación estándard de la población es calculada
                  na.rm = T # Si es verdadero (T), entonces valores NA son ignorados
                  )

¿Es lo esperado?

plot(tmean_sd)

Tips!

Recuerde que en RStudio puede iluminar la función (ej., stdev) y oprimir F1 para obtener ayuda, es decir información sobre los argumentos que se encuentran dentro de una función.

  • Calcular la temperatura media utilizando app()

También podemos calcular la temperatura media del año, utilizando la función app() del paquete terra.

tmean_mean <- app(x = tmean,  # SpatRaster
                  fun = mean) # función
plot(tmean_mean, main="Temperatura media anual (°C)")

Métodos zonales y globales

Los métodos zonales y globales hacen referencia a hacer cálculos del area total del ráster utilizando la función global() o zonas de un ráster utilizando la función zonal() (estas zonas se pueden definir con otro capa). Por ejemplo, la figura muestra cómo podemos calcular la sumatoria de una zona en el ráster.

Veamos un ejemplo de cómo podemos utilizar la base de datos global de la huella humana (Human footprint, en inglés) para calcular el promedio de huella humana para cada país del mundo.

La huella humana reune la mayor parte de las acciones humanas con el potencial de deteriorar la naturaleza (ver el trabajo de Venter et al 2016 para mayores detalles).

La escala es 0 = muy baja, hasta 50 = muy alta huella humana.

Note que en nuestros códigos vamos a utilizar el simbolo %>%, conocido como forward-pipe operator en inglés, que envia el resultado de una expresión a la otra función dentro de un mismo objeto. %>% fue originalmente desarrollado paquete magrittr.

En este caso cargamos la base de datos “footprint” y después la agregamos 10 veces utilizando la función aggregate(). Tenga en cuenta que la función aggregate() se encuentra en varios paquetes, entonces es buena práctica especificar que estámos utilizando la función del paquete terra, asi: terra:aggregate().

¿Por qué agregamos el ráster huella humana? Vamos a convertir la huella humana en un raster de resolución más gruesa para acelerar el procesamiento de nuestros siguientes pasos (es decir menos tiempo de procesamiento). Más adelante veremos unos ejemplos del tiempo de procesamiento en diferentes resoluciones.

huella <- geodata::footprint(year = 2009,
                             path ="C:/Talleres_Ecodiversa_R/Intro_espacial")%>% 
          terra::aggregate(fact=10, fun=mean)
## 
|---------|---------|---------|---------|
=========================================
                                          
huella
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             45.21

Ahora carguemos los países del mundo.

paises <- geodata::world(resolution = 5,
                         level = 0,
                        path = "C:/Talleres_Ecodiversa_R/Intro_espacial")
paises
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 231, 2  (geometries, attributes)
##  extent      : -180, 180, -90, 83.65625  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
##  names       : GID_0      NAME_0
##  type        : <chr>       <chr>
##  values      :   ABW       Aruba
##                  AFG Afghanistan
##                  AGO      Angola

Es una base de datos que contiene información vectorial. Más adelante dedicaremos una sección a este formato, por ahora transformemos la información vectorial a ráster, utilizando la función rasterize() del paquete terra.

Note que vamos a utilizar el ráster huella como plantilla para rasterize() los polígonos paises; es decir obtendremos un ráster de las mismas dimensions:, resolution:, extent y coord. ref.:. Para los valores del ráster utilizaremos el nombre de la columna en el objeto vectorial que contiene el nombre de los países, es decir NAME_0.

paises_r <- terra::rasterize(paises, huella, field="NAME_0")
paises_r
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :      NAME_0 
## min value   : Afghanistan 
## max value   :       Åland
plot(paises_r)

Ahora si, apliquemos la función zonal() del paquete terra para calcular el promedio de la huella humana en cada país del mundo.

huella_paises <- terra::zonal(huella,
                              paises_r,
                              as.raster=T, # si es F obtiene solo los datos.
                              fun = mean
                              )
plot(huella_paises)

Seleccionado un área geográfica dentro de un ráster

Es muy común que usted tenga un ráster con una extensión global, pero usted quiera enfocar su análisis en una región específica. Para esto vamos a utilizar la función crop() del paquete terra.

Por ejemplo, seleccionemos los datos de la huella humana solo para Colombia.

Las categorias del objeto paises_r pueden ser inspeccionadas con la función levels() del paquete terra. Generalmente son dos columnas, la primera es un numero entero que identifica los valores de la celda y la segunda es una etiqueta que indica el nombre de la categoria. Si revisa la tabla podrá ver que Colombia es el valor 43.

knitr::kable(head(terra::levels(paises_r)))%>%
  kableExtra::kable_paper("hover", full_width = F) %>%
  kableExtra::scroll_box(width = "100%", height = "200px")
value NAME_0
0 Afghanistan
1 Akrotiri and Dhekelia
2 Albania
3 Algeria
4 American Samoa
5 Andorra
6 Angola
7 Antarctica
8 Antigua and Barbuda
9 Argentina
10 Armenia
11 Aruba
12 Australia
13 Austria
14 Azerbaijan
15 Bahamas
16 Bahrain
17 Bangladesh
18 Barbados
19 Belarus
20 Belgium
21 Belize
22 Benin
23 Bhutan
24 Bolivia
25 Bonaire, Saint Eustatius and Saba
26 Bosnia and Herzegovina
27 Botswana
28 Brazil
29 Brunei
30 Bulgaria
31 Burkina Faso
32 Burundi
33 Cabo Verde
34 Cambodia
35 Cameroon
36 Canada
37 Caspian Sea
38 Cayman Islands
39 Central African Republic
40 Chad
41 Chile
42 China
43 Colombia
44 Comoros
45 Congo
46 Cook Islands
47 Costa Rica
48 Croatia
49 Cuba
50 Curaçao
51 Cyprus
52 Czech Republic
53 Côte d’Ivoire
54 Democratic Republic of the Congo
55 Denmark
56 Djibouti
57 Dominica
58 Dominican Republic
59 East Timor
60 Ecuador
61 Egypt
62 El Salvador
63 Equatorial Guinea
64 Eritrea
65 Estonia
66 Eswatini
67 Ethiopia
68 Falkland Islands
69 Faroe Islands
70 Fiji
71 Finland
72 France
73 French Guiana
74 French Polynesia
75 French Southern Territories
76 Gabon
77 Gambia
78 Georgia
79 Germany
80 Ghana
81 Greece
82 Greenland
83 Grenada
84 Guadeloupe
85 Guam
86 Guatemala
87 Guinea
88 Guinea-Bissau
89 Guyana
90 Haiti
91 Heard Island and McDonald Islands
92 Honduras
93 Hong Kong
94 Hungary
95 Iceland
96 India
97 Indonesia
98 Iran
99 Iraq
100 Ireland
101 Isle of Man
102 Israel
103 Italy
104 Jamaica
105 Japan
106 Jersey
107 Jordan
108 Kazakhstan
109 Kenya
110 Kiribati
111 Kosovo
112 Kuwait
113 Kyrgyzstan
114 Laos
115 Latvia
116 Lebanon
117 Lesotho
118 Liberia
119 Libya
120 Liechtenstein
121 Lithuania
122 Luxembourg
123 Macao
124 Macedonia
125 Madagascar
126 Malawi
127 Malaysia
128 Mali
129 Malta
130 Martinique
131 Mauritania
132 Mauritius
133 Mayotte
134 Mexico
135 Micronesia
136 Moldova
137 Mongolia
138 Montenegro
139 Montserrat
140 Morocco
141 Mozambique
142 Myanmar
143 Namibia
144 Nepal
145 Netherlands
146 New Caledonia
147 New Zealand
148 Nicaragua
149 Niger
150 Nigeria
151 Niue
152 North Korea
153 Northern Cyprus
154 Northern Mariana Islands
155 Norway
156 Oman
157 Pakistan
158 Palau
159 Palestine
160 Panama
161 Papua New Guinea
162 Paraguay
163 Peru
164 Philippines
165 Poland
166 Portugal
167 Puerto Rico
168 Qatar
169 Romania
170 Russia
171 Rwanda
172 Réunion
173 Saint Helena
174 Saint Kitts and Nevis
175 Saint Lucia
176 Saint Pierre and Miquelon
177 Saint Vincent and the Grenadines
178 Saint-Martin
179 Samoa
180 Saudi Arabia
181 Senegal
182 Serbia
183 Seychelles
184 Sierra Leone
185 Singapore
186 Sint Maarten
187 Slovakia
188 Slovenia
189 Solomon Islands
190 Somalia
191 South Africa
192 South Georgia and the South Sandwich Islands
193 South Korea
194 South Sudan
195 Spain
196 Sri Lanka
197 Sudan
198 Suriname
199 Svalbard and Jan Mayen
200 Sweden
201 Switzerland
202 Syria
203 São Tomé and Príncipe
204 Taiwan
205 Tajikistan
206 Tanzania
207 Thailand
208 Togo
209 Tonga
210 Trinidad and Tobago
211 Tunisia
212 Turkey
213 Turkmenistan
214 Turks and Caicos Islands
215 Uganda
216 Ukraine
217 United Arab Emirates
218 United Kingdom
219 United States
220 Uruguay
221 Uzbekistan
222 Vanuatu
223 Venezuela
224 Vietnam
225 Virgin Islands, U.S.
226 Western Sahara
227 Yemen
228 Zambia
229 Zimbabwe
230 Åland


Lo que vamos hacer es copiar el objeto paises_r y llamarlo Colombia, luego vamos a reclasificar todos los valores diferentes a 43 con NA (sin valores), es decir escoger Colombia.

Colombia <- paises_r
Colombia[Colombia != 43]<- NA
Colombia
## class       : SpatRaster 
## dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia

Visualicemos el mapa de Colombia. Hemos seleccionado lo que queremos, pero el extent: sigue siendo el mismo, es decir todo el globo entero, por eso vemos el país tan pequeño.

plot(Colombia)

Podemos utilizar la función trim() para remover las filas y columnas con valores NA (alrededor de los valores del mapa de Colombia), y como consecuencia se cambia el extent: del nuevo objeto.

Colombia_t <- terra::trim(Colombia)
Colombia_t
## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : NAME_0 
## name        :   NAME_0 
## min value   : Colombia 
## max value   : Colombia
plot(Colombia_t)

Ahora si utilicemos la función crop() para extraer los datos de temperatura media en el extent del objeto Colombia_t.

Colombia_huella <- terra::crop(huella, Colombia_t)
Colombia_huella
## class       : SpatRaster 
## dimensions  : 200, 146, 1  (nrow, ncol, nlyr)
## resolution  : 0.08333333, 0.08333333  (x, y)
## extent      : -79, -66.83333, -4.25, 12.41667  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : wildareas-v3-2009-human-footprint 
## min value   :                              0.00 
## max value   :                             37.86

Visualicemos el nuevo objeto. Como podemos ver la función crop() extrae lo valores de la huella humana con base en el extent: del objeto Colombia_t.

plot(Colombia_huella)

Para tener solamente los valores de la huella humana dentro de los limites de Colombia_t debemos utilizar la función mask(). Veamos:

Colombia_huella_mask <- terra::mask(Colombia_huella, Colombia_t)
plot(Colombia_huella_mask)

Ahora bien, por qué no utilizar la opción mask() desde un principio? Hay varias razones.

  • las dos raster deben ser del mismo extent:. Al intentar implementar la función mask() directamente al mapa global huella obtenemos el Error. Es decir el mapa huella (es global) y el mapa Colombia_t tienen diferente extent:.
Colombia_huella_2 <- terra::mask(huella, Colombia_t)
## Error: [mask] extents do not match
  • la segunda razón es que es más eficiente para R implementar la función crop(), porque selecciona filas y columnas del área definida por el extent: que queremos. Al contrario, mask() require hacer más cálculos para seleccionar los pixeles que estan en la figura compleja del mapa de Colombia.

Tips!

Por eso, es importante implementar primero crop() y luego mask() para acelerar el proceso.

Seleccionando valores en un ráster

En ocasiones reqerimos utilizar solo unos valores determinados dentro del ráster.

Por ejemplo, queremos obtener un subset de las temperaturas por encima de 25 grados. Para lo cual podemos reclasificar los valores menores a 25 con NA (no values).

tmean_caliente <- tmean_mean
tmean_caliente[tmean_caliente < 25] <- NA
tmean_caliente
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :    mean 
## min value   : 25.0000 
## max value   : 30.9879
plot(tmean_caliente)

Pero tambien podemos necesitar un rango , por ejemplo entre 15 y 20 grados. Para eso utilizamos el operador | que indica “o” (los valores menores a 15 “o” valores mayores a 20, serán iguales a NA).

tmean_templado <- tmean_mean
tmean_templado[tmean_templado < 15 | tmean_templado >20] <- NA
plot(tmean_templado)

Tambien podemos extraer lo valores en una pila de rásters. Por ejemplo

tmean_frio <- tmean
tmean_frio[tmean_frio > 10] <- NA
tmean_frio
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## names       : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ... 
## min values  :     -45.884,       -44.8,   -57.92575,    -64.1925,    -64.8115,   -64.35825, ... 
## max values  :      10.000,        10.0,    10.00000,     10.0000,     10.0000,    10.00000, ...

Seleccionando valores en un ráster con otro ráster

Vamos a extraer los valores de precipitación global en areas con temperaturas mayores a 20°C (tmean_caliente).

Primero, descarguemos la precipitación utilizando el paquete geogata:

prec <- geodata::worldclim_global(
          var="prec",
          res=10,
          path = "C:/Talleres_Ecodiversa_R/Intro_espacial"
                                 )
prec
## class       : SpatRaster 
## dimensions  : 1080, 2160, 12  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : wc2.1_10m_prec_01.tif  
##               wc2.1_10m_prec_02.tif  
##               wc2.1_10m_prec_03.tif  
##               ... and 9 more source(s)
## names       : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :         908,         793,         720,        1004,        2068,        2210, ...

Ahora calculemos la precipitación media anual:

prec_mean <- app(prec, mean)
prec_mean
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833
plot(prec_mean)

Finalmente, utilicemos la función mask() para extraer los valores de precipitación en areas con temperaturas mayores a 20°C.

prec_caliente <- terra::mask(prec_mean, tmean_caliente)
prec_caliente
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     mean 
## min value   :   0.0000 
## max value   : 932.5833

Visualizemos los dos rásters:

# Creemos dos columnas para los plots
par(mfrow=c(1,2))

# Temperature > 20 °C
plot(tmean_mean*0, legend=F) # Fondo gris para mostrar los continentes 
t_breaks <- seq(minmax(tmean_caliente)[1], minmax(tmean_caliente)[2], 1) # escala de la temperatura
terra::plot(tmean_caliente, add=T, type="interval", plg=list(x="bottom", cex=1.2, title="°C"), breaks = t_breaks)
title("Temperatura media > 20 °C", line = -2, cex.main=1.6)

# Precipitación en áreas con temp  > 20 °C
p_breaks <- seq(minmax(prec_caliente)[1], minmax(prec_caliente)[2], 150)
plot(tmean_mean*0, legend=F)
terra::plot(prec_caliente,  add=T, type="interval", plg=list(x="bottom", cex=1.2, title="mm"), col=rev(topo.colors(15)), breaks=p_breaks)
title("Precipitación en áreas con temp  > 20 °C", line = -2, cex.main=1.6)

Extraer los valores de un ráster

En algún momento de sus análisis necesitará los valores de los rásters para hacer algunos cálculos adicionales o preparar gráficas. El paquete terra ofrece la función values() para extraer lo valores del ráster.

Veamos un ejemplo con la temperatura media y la precipitación global.

Visualicemos una gráfica que muestre los valores pareados de temperatura media y precipitación.

Interesante ver que las precipitaciones medias anuales en áreas frías son menores que áreas más calientes.

plot(values(tmean_mean), values(prec_mean))

Ahora tomemos solamente los datos de precipitación para Colombia.

Colombia_t_prec <- crop(prec_mean,  # precipitación media
                        aggregate(Colombia_t, 2) # agreguemos colombia 2 veces para que sea la misma resolución de la precipitación.
                        )%>%
                  resample(aggregate(Colombia_t, 2))%>% # asegurarnos que las raster sean iguales en extent, resolution, etc.
                  mask(aggregate(Colombia_t, 2))       # Seleccionar solamente el borde administrativo de Colombia

plot(Colombia_t_prec)

El Pacífico siempre lluvioso!

Seleccionemos también la temperatura media para Colombia.

Colombia_t_tmean <- crop(tmean_mean, aggregate(Colombia_t, 2))%>%
                    resample(aggregate(Colombia_t, 2))%>%
                    mask(aggregate(Colombia_t, 2))
plot(Colombia_t_tmean)

Ahora visualicemos los valores pareados de temperatura media y precipitación para Colombia con relación al global.

Definitivamente tropical!

  plot(values(tmean_mean), values(prec_mean), xlim=c(-55, 30), ylim=c(0, 1000),
       ylab="",yaxt="n", xlab="",xaxt="n", pch=21, col="black")

par(new=T)
plot(values(Colombia_t_tmean), values(Colombia_t_prec), col="red", xlim=c(-55, 30), ylim=c(0, 1000), pch=22)

legend("topleft", legend=c("Colombia", "Global"), col=c("red", "black"), pch=c(22,21), cex=0.8)

Proyectando rásters

Hasta ahora hemos utilizado coordenadas geográficas para estos análisis. Sin embargo, los análisis más locales o regionales requieren de ráster proyectadas. Es decir los datos proyectados en una superficie plana. En especial si queremos hacer cálculos de áreas para informar sobre aspectos del área de estudio (ej., área de bosques fragmentados, etc).

Hay diferentes proyecciones cartográficas, pero la que más vamos a utilizar en el análisis espacial es la que mantienen el sentido del área. Cada país y región tiene proyecciones apropiadas que preservan el area, entonces hay que hacer el ejercicio de buscar la mejor proyección para el área de estudio.

Por fortuna hay muchos sitios para buscar y spatialrefence y epsg.io.

Por ejemplo, si buscamos “Colombia” en “spatialreference” nos aparecerá algo asi:

En nuestros códigos podemos utilizar cualquiera de los dos formatos ESPG (es el que aparece en el primer pantallazo de ‘spatialreference’) o Proj4.

Por ejemplo si utilizamos EPSG:3116 el equivalente Proj4 es: +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Veamos cómo utilizamos esto en R.

  • Carguemos la cobertura vegetal (solamente la clase árboles (trees)),
  • seleccionemos los datos para Colombia y
  • utilicemos la proyección EPSG:3116 para reproyectar el ráster.

La función project() de terra incluye varios argumentos:

  • y = hace referencia a “coordinate reference system, (por sus siglas en inglés), en este caso EPSG:3116.
  • method = “bilinear” es util para variables continuas y “near” para variables categóricas. Vamos a utilizar “bilinear porque es la proporción de bosque es cada celda (es decir que es una variable continua).
  • res = hace referencia al tamaño de celda (en metros) que queremos al final. Es decir 1Km por 1Km. Celdas cuadradas, lo que facilitará y mejorará la exactitud (accuracy, en inglés) de los cálculos.
arboles <- geodata::landcover(var= 'trees',
                    path = "C:/Talleres_Ecodiversa_R/Intro_espacial"                            )%>%
  terra::crop(disagg(Colombia_t,2))%>%
  terra::resample(terra::disagg(Colombia_t,2))%>%
  terra::mask(disagg(Colombia_t,2))%>%
  terra::project(y="EPSG:3116", method="bilinear", res = 1000)
  
arboles
## class       : SpatRaster 
## dimensions  : 1857, 1358, 1  (nrow, ncol, nlyr)
## resolution  : 1000, 1000  (x, y)
## extent      : 451351.3, 1809351, 19203.33, 1876203  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Colombia Bogota zone (EPSG:3116) 
## source(s)   : memory
## name        : trees 
## min value   :     0 
## max value   :     1
plot(arboles)

Utilicemos la función cellSize() de terra para calcular el área de cada celda de este ráster. Casi todas las celdas son cercanas a 1Km2, pero hay algunas muy pequeñas variaciones, en especial en el oriente de Colombia, porque esta proyección esta centrada más en la mitad de Colombia.

Continuemos con nuestro análisis, asumiendo que las celdas son de igual tamaño. Sin embargo, tenga en cuenta reportar la exactitud en sus resultados.

arboles_cellsize <- cellSize(arboles, unit="km")
plot(arboles_cellsize)

Reclasificando rásters

La reclasificación de rásters es de mucha utilidad en el análisis espacial. Por ejemplo si usted quiere reclasificar la cobertura de arboles en Colombia en cuatro categorias.

Para esto debemos crear una matriz que contenga:

  • En las columnas el valor inicial, el valor final del raster y que valor quermos al final

  • Las filas son determinadas con base en los rangos que usted defina.

Por ejemplo, definamos tres nuevas categorías de acuerdo a la proporción de árboles y reclasifiquemolas con 1, 2 y 3. En la primera linea se dice, de 0 a 0.25 reclasifiquelo a 1.

m <- c(0, 0.25, 1,
       0.25, 0.5, 2,
       0.5, 1, 3)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
rclmat
##      [,1] [,2] [,3]
## [1,] 0.00 0.25    1
## [2,] 0.25 0.50    2
## [3,] 0.50 1.00    3

Ahora apliquemos esta matriz al ráster.

rc <- classify(arboles,
               rclmat)
plot(rc, legend = FALSE,
     col = c("dark green",  "red", "yellow"), axes = FALSE)
legend("topright",
       legend = c("Alta", "Media",  "Baja"),
       fill = c("dark green", "yellow",  "red"),
       border = FALSE,
       bty = "n") # turn off legend border

Podríamos tener unas datos del área de cada categoría de la siguiente forma:

terra::freq(rc)
##   layer value  count
## 1     1     0     91
## 2     1     1 213491
## 3     1     2 164580
## 4     1     3 761372

Tiempo de procesamiento

pixel <- c(1, 10, 100, 1000)

pixel_size <- list()
for(i in 1:length(pixel)){
  
pixel_size[[i]] <- rast(ncol=pixel[i], nrow=pixel[i])
  values(pixel_size[[i]]) <- 1:ncell(pixel_size[[i]])

}
elevation <- geodata::elevation_global(res = 0.5,
                                       path = "C:/Talleres_Ecodiversa_R/Intro_espacial")
elevation
## class       : SpatRaster 
## dimensions  : 21600, 43200, 1  (nrow, ncol, nlyr)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : wc2.1_30s_elev.tif 
## name        : wc2.1_30s_elev 
## min value   :           -415 
## max value   :           8424
plot(elevation)

elevation_2000m <-  elevation
system.time(elevation_2000m[elevation_2000m < 2000] <- NA)
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
##    user  system elapsed 
##   74.25    4.42   78.68
plot(elevation_2000m)

Cargando y guardando rasters

Formatos de datos en ráster

Existen varios tipos o formatos de datos en ráster, como por ejemplo GeoTIFF, Esr grid, JPEG 2000, ERDAS IMAGINE (IMG), GeoPackage MBTiles, ASCII, entre otros.

Puede encontrar una descripció detallada de este tipo de formatos de ráster en este sitio.

La mayoría de formatos están disponibles en R tanto para leer el ráster (raster o terra) como para guardar en raster o en terra.

8. Vectorial

11. Análisis espaciales aplicados a la conservación de la biodiversidad.

Análisis globales vs. locales (el problema de escala)

Cambios en la cobertura vegetal

Cambios en la dsitribuión de las especies

12. Informacion adicional

Para incluir…………… Sin embargo, existen muchas bases de datos que pueden ser importados a un formato vectorial .Por ejemplo, registros que tengan coordenadas geográficas en una hoja de cálculo o texto podrían ser importados y crear un objeto vectorial (ej., puntos).